Conserving GW scheme for nonequilibrium quantum transport in molecular contacts 
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We give a detailed presentation of our recent scheme to include correlation effects in molecular transport 
calculations using the non-equilibrium Keldysh formalism. The scheme is general and can be used with any 
quasiparticle self-energy, but for practical reasons we mainly specialize to the so-called GW self-energy, widely 
used to describe the quasiparticle band structures and spectroscopic properties of extended and low-dimensional 
systems. We restrict the GW self-energy to the central region, and describe the leads by density functional the- 
ory (DFT). A minimal basis of maximally localized Wannier functions is applied both in the central GW region 
and the leads. The importance of using a conserving, i.e. fully self-consistent, GW self-energy is demonstrated 
both analytically and by numerical examples. We introduce an effective spin-dependent interaction which auto- 
matically reduces self-interaction errors to all orders in the interaction. The scheme is applied to the Anderson 
model in- and out of equilibrium. In equilibrium at zero temperature we find that GW describes the Kondo 
resonance fairly well for intermediate interaction strengths. Out of equilibrium we demonstrate that the one- 
shot GqWo approximation can produce severe errors, in particular at high bias. Finally, we consider a benzene 
molecule between featureless leads. It is found that the molecule's HOMO-LUMO gap as calculated in GW 
is significantly reduced as the coupling to the leads is increased, reflecting the more efficient screening in the 
strongly coupled junction. For the IV characteristics of the junction we find that HF and GqWoIGhf] yield 
results closer to GW than does DFT and GoVKo[Gdft]. This is explained in terms of self-interaction effects and 
life-time reduction due to electron-electron interactions. 

PACS numbers: 72.10.-d,71.10.-w,73.63.-b 



I. INTRODUCTION 

Since the first measurements of electron transport through 
single molecules were reported in the late nineties the the- 
oretical interest for quantum transport in nano-scale systems 
has been rapidly growing. An important driving force behind 
the scientific developments is the potential use of molecular 
devices in electronics and sensor applications. On the other 
hand it is clear that a successful introduction of these tech- 
nologies is heavily dependent on the availability of theoreti- 
cal and numerical tools for the accurate description of such 
molecular devices. 

So far, the combination of density functional theory (DFT) 
and non-equilibrium Green's functions (NEGF) has been 
the most popular method for modeling nano-scale conduc- 
tivity'*'^'^-^. For strongly coupled systems such as metallic 
point contacts, monatomic chains, and contacts with small 
chemisorbed molecules, this combination has been remark- 
ably successful^'^'"', but in the opposite limit of weakly cou- 
pled systems where the conductance is much smaller than the 
conductance quantum, Gq = 2e^/h, the NEGF-DFT method 
has been found to overestimate the conductance relative to ex- 
periments' Part of this discrepancy might result from 
the use of inappropriate exchange-correlation (xc) function- 
als'"*. However, it is important to remember that the applica- 
tion of ground state DFT to non-equilibrium transport cannot 
be rigorously justified - even with the exact xc-functional. In 
particular, a breakdown of the effective single-particle DFT 
description is expected when correlation effects are important 
or when the system is driven out of equilibrium. 



Over the years several different schemes have been pro- 
posed as alternatives to NEGF-DFT. Historically, the first 
DFT based transport methods used an equivalent formula- 
tion in terms of scattering states rather than Green's func- 
jjQjjgi5ji6J2 A more recent approach (still within DFT), solves 
a master equation for the density matrix of an electron system 
exposed to a constant electric field and coupled to a damping 
heat bath of auxiliary phononsi^. 

A few attempts have been made to calculate the current in 
the presence of electronic correlations. In one approach the 
density matrix is obtained from a many-body wave function 
and the non-equilibrium boundary conditions are invoked by 
fixing the occupation numbers of left- and right going states'-. 
Exact diagonalization within the molecular subspace has been 
combined with rate equations to calculate tunneling currents 
to first order in the lead-molecule coupling strength?^. The lin- 
ear response conductance of jellium quantum point contacts 
has been addressed on the basis of the Kubo formula^''^^. Al- 
though this method is restricted to the low bias regime, it has 
the advantage over the NEGF method that interactions outside 
the device region can be naturally included. The time depen- 
dent version of density functional theory has also been used as 
framework for quantum transport^^'^'*-^^. This scheme is par- 
ticularly useful for simulating transients and high frequency 
ac -responses. Within the NEGF formalism the many-body 
GW approximation has been used to address coiTelated trans- 
port both under equilibrium^- and non-equilibriuro21 condi- 
tions. 

Within the framework of many-body perturbation theory 
(MBPT) electronic correlations are described by a self-energy 
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which in practice must be obtained according to some approx- 
imate scheme, e.g. by summing a restricted set of Feynman 
diagrams. The important question then arises whether the 
quantities calculated from the resulting Green's function will 
obey the simple conservation laws. In the context of quan- 
tum transport the continuity equation, which ensures charge 
conservation, is obviously of special interest. An elegant way 
of invoking the conservation laws is to write the self-energy 
as the functional derivative of a so-called ^-functional, i.e. 
S[G] = 5<^[G]/5G. Since the self-energy in this way be- 
comes dependent on the Green's function (GF), it must be 
determined self-consistently in conjunction with the Dyson 
equation. 

Due to the large computational demands connected with 
the self-consistent solution of the Dyson equation, practical 
GW band structure calculations usually evaluates the self- 
energy at some approximate non-interacting Go. This non 
self-consistent scheme does not constitute a conserving ap- 
proximation. While this might not be important for the cal- 
culated spectrum, self-consistency has been demonstrated to 
be fundamental for out-of-equilibrium transport^^. In addi- 
tion to its conserving nature, another nice feature of the self- 
consistent approach is that it leads to a unique GF and thus re- 
moves the Go -dependence inherent in the non self-consistent 
approach. 

A reliable description of electron transport through a 
molecular junction requires first of all a reliable description of 
the internal electronic structure of the molecule itself, i.e. its 
electron addition and removal energies. The GW approxima- 
tion has been widely and successfully used to calculate such 
quasiparticle excitations in both semi conductors, insulators 
and molecules^2i^2i2L22i22, and on this basis it seems natural to 
extend its use to transport calculations. 

There are two main obstacles related to the extension of the 
GW method to charge transport. First, the conventional appli- 
cation of the GW method has been to ground state problems 
whereas transport is an inherent non-equilibrium problem. 
Secondly, it is not obvious how to treat electron-electron in- 
teractions in the leads within the NEGF formalism. In Ref.I^^ 
we proposed to overcome these problems by extending the 
GW self-energy to the Keldysh contour and restricting it to 
a finite central region where correlation effects are expected 
to be most important. In the present paper we provide an ex- 
tended presentation of these ideas. 

When a molecule is brought into contact with electrodes 
a number of physical mechanisms will affect its electronic 
structure. Some of these mechanisms are single-particle in 
nature and are already well described at the DFT Kohn- 
Sham level. But there are also important many-body ef- 
fects which require a dynamical treatment of the electronic 
interactions. One example is the renormalization of the 
HOMO-LUMO gap induced by the image charges formed 
in the electrodes when an electron is added to or removed 
from the moleculei^^i^ Another example is the Kondo effect 
which results from correlations between a localized spin on 
the molecule and delocalized electrons in the electrodes^^'^*. 
Third, as we will show here, the coupling to (non-interacting) 
electrodes enhances the screening on the molecule leading to 



acharacteristic reduction of the HOMO-LUMO gap as func- 
tion of the electrode-molecule coupling strength. 

In this paper, we focus on improving the description of 
quantum transport in molecular junctions by improving the 
description of the internal electronic structure of the molecule 
while preserving a non-perturbative treatment of the coupling 
to leads. We do this within the NEGF formalism by us- 
ing a self-consistent GW self-energy to include xc effects 
within the molecular subspace which in turn is coupled to non- 
interacting leads. The rationale behind this division is that the 
transport properties to a large extent are determined by the 
narrowest part of the conductor, i.e. the molecule, while the 
leads mainly serve as particle reservoirs. Strictly speaking this 
is correct only when a sufficiently large part of the leads is in- 
cluded in the GW region. If the central region is too small, 
spurious back-scattering at the interface between the GW and 
mean-field regions might affect the calculated conductance. 
Furthermore, the dynamical formation of image charges in the 
electrodes requires that part of the electrodes are included in 
the GW region. In the present work we do, however, not at- 
tempt to address this latter effect. 

The paper is organized as follows. In Sec. |II]we introduce 
the model used to describe the transport problem and review 
the basic elements of the Keldysh Green's function formalism. 
In Sec.|III]we introduce an effective interaction, discuss the 
problem of self-interaction correction in diagrammatic expan- 
sions, and derive the non-equilibrium GW equations for an 
interacting region coupled to non-interacting leads. In Sec. II VI 
we introduce the current formula and show that charge con- 
servation is fulfilled within the NEGF formalism for $ deriv- 
able self-energies - also when incomplete basis sets are used. 
The practical implementation of the GW transport scheme us- 
ing a Wannier function basis obtained from DFT is described 
in Sec. |V] In Sees. [VT] and IVIII we present results for the 
non-equilibrium transport properties of the Anderson impu- 
rity model and benzene molecule between jellium leads, re- 
spectively. We conclude in Sec. I Villi 



II. GENERAL FORMALISM 

In this section we review the elements of the Keldysh 
Green's function formalism necessary to deal with the non- 
equilibrium transport problem. To limit the technical details 
we specialize to the case of orthogonal basis sets and refer to 
Ref . I37I for a generalization to the non-orthogonal case. 



A. Model 

We consider a quantum conductor consisting of a central re- 
gion (G) connected to left {L) and right (i?) leads. For times 
t < to the three regions are decoupled from each other, each 
being in thermal equilibrium with a common temperature, T, 
and chemical potentials fiL,IJ-c, ™d /i^, respectively. At 
t ^ to the coupling between the three subsystems is switched 
on and a current starts to flow as the electrode with higher 
chemical potential discharges through the central region into 
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the lead with lower chemical potential. Our aim is to calculate 
the steady state current which arise after the transient has died 
out. 

We denote by an orthonormal set of single-particle 
orbitals, and by H the Hilbert space spanned by {(f>i}. The 
orbitals (pi are assumed to be locaUzed such that H can be de- 
composed into a sum of orthogonal subspaces corresponding 
to the division of the system into leads and central region, i.e. 
H = Hl + 'He + H-R. We will use the notation i e a to 
indicate that (pi £ TL^ for some a G {i, C, R}. 

The non-interacting part of the Hamiltonian of the con- 
nected system is written 



E E 

'.je o-=Ti 

L.CR 



(1) 



where i,j run over all basis states of the system. For a, f3 € 
{L, C, R}, the operator hap is obtained by restricting i to re- 
gion a and j to region (3 in Eq. ([T]i. Occasionally we shall 
write ha instead of haa- We assume that there is no direct 
coupling between the two leads, i.e. h^n — hjiL = (this 
condition can always be fulfilled by increasing the size of the 
central region since the basis functions are localized). We in- 
troduce a special notation for the "diagonal" of h. 



ho = hLL + hoc + h 



RR- 



(2) 



It is instructive to note that h^ does not describe the three re- 
gions in isolation from each other, but rather the contacted 
system without inter-region hopping. We allow for interac- 
tions between electrons inside the central region. The most 
general form of such a two-body interaction is. 



V 



(3) 



replacing L by R. For qc and Zc we must add V to account 
for correlations in the initial state of the central region. The 
initial state of the whole system is then given by 



Q = QlQcQR, 



(7) 



If V is not included in qc we obtain the uncorrelated (non- 
interacting) initial state Qni ■ We note that the order of the den- 
sity matrices in Eq. (|7| plays no role since they all commute 
due to the orthogonality of the system {(pi}. Because (ho) 
describes the contacted system without inter-region hopping, 
g {Qni) does not describe the three regions in physical isola- 
tion. In other words the three regions are only decoupled at 
the dynamic level for times t < to. 



Left lead 
(L) 



Central region 
(C) 



Right lead 

(R) 



4^ 
4^ 



FIG. 1: Before the coupling between the three regions is established, 
the three subsystmes are in equilibrium with chemical potentials hl, 
fic, and respectively. 



The full Hamiltonian describing the system at time t can then 
be written 



H{t) 



Ho = ho + V for t < to 
H = h + V for t > to 



(4) 



Notice, that we use small letters for non-interacting quantities 
while the subscript refers to uncoupled quantities. The spe- 
cific form of the matrix elements hij and Vij^ki defining the 
Hamiltonian are considered in Sec. |V] 

Having defined the Hamiltonian we now consider the intial 
state of the system, i.e. the state at times t < to. For such 
times the three subsystems are each in thermal equilibrium 
and thus characterized by their equilibrium density matrices. 
For the left lead we have 
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with 



0L 



Zr 



Tr[exp(-/3(/iL-AiL7VL))]. 



Here j3 is the inverse temperature and = Tla 
is the number operator of lead L. qr and Zr are obtained by 



(5) 



(6) 

^irr^it^ 



B. The Contour-ordered Green's function 

In this section we introduce the contour-ordered GF which 
is the central object for the many-body perturbation theory in 
non-equilibrium systems. For more detailed accounts of the 
NEGF theory we refer to Refs. |38ll39i 

The contour-ordered GF relevant for the model introduced 
in the previous section, is defined by 

G..,j.' (r, r') = -lTy{QT[cH.^a[T)c^H,,.' (^')]}- (8) 

Here r and r' are points on the Keldysh contour, C, which 
runs along the real time axis from to to oo and back to to, and 
T is the time-ordering operator on the contour. The creation 
and annihilation operators are taken in the Heisenberg picture 
with respect to the full Hamiltonian in Eq. (|4]i. We do not 
consider spin-flip processes and thus suppress the spin indices 
in the following. 

In order to obtain an expansion of Gij (r, r') in powers of 
V, we switch to the interaction picture where we have 

G,,(r,r') = ~iTr{QT[e-^ic^^^-^^)cr,,,{T)cl ,(t')]}. (9) 
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By extending C into the complex plane by a vertical branch 
running from to — if], we can replace g by the uncorre- 
cted Qm^"^- Neglecting the vertical branch then corresponds 
to neglecting correlations in the central region's initial state. 
While it must be expected that the presence of initial corre- 
lations will influence the transient behavior of the current, it 
seems plausible that they will be washed out over time such 
that the steady state current will not depend on gc- Further- 
more, in the special case of equilibrium (/i^ — jic — fJ.R) 
and zero temperature, the Gellman-Low theorem ensures that 
the correlations are correctly introduced when starting from 
the uncorrected initial state at to = — 00.**" In practice the 
neglect of initial correlations is a major simplification which 
allows us to work entirely on the real axis avoiding any refer- 
ence to the imaginary time. For these reasons we shall adopt 
this approximation and neglect initial correlations in the rest 
of this paper. 

Eq. (|9]l with g replaced by gni constitute the starting point 
for a systematic series expansion of Gij in powers of V and 
the free propagator, 

5„(r,r') = -iTr{g„,T[chAr)cl^ir')]}, (10) 

which describes the non-interacting electrons in the coupled 
system. The diagrammatic expansion leads to the identifica- 
tion of a self-energy, E, which relates the interacting GF to 
the non-interacting one through Dyson's equation 

G{T,T')=g{T,T')+ f dridT2<7(r,ri)E(Ti,r2)G(r2,r'), 
Jc 

(11) 

(matrix multiplication is implied). As we will see in 
Sec. IIV Al only the Green's function of the central region is 
needed for the calculation of the current, and we can therefore 
focus on the central-region submatrix of G. Due to the struc- 
ture of V, the self-energy matrix, Ey , will be non-zero only 
when both i,j g G, and for this reason G subscripts can be 
added to all matrices in Eq. ( fTTT ). Having observed this we will 
nevertheless write S instead of for notational simplicity. 

The free propagator, gc{T,T'), which is still a non- 
equilibrium GF, satisfies the following Dyson equation 

9c{r, t') = 5o,c(t, t') + J dridr25o,c(T, n) 

[I]L(ri,T2) + I]fl(ri,T2)]5c(T-2,r'), (12) 

where go is the equilibrium GF defined by gni and ho- The 
coupling self-energy due to lead a ^ L,R is given by 

Sa(T,T') = hcagQ,aiT,T')haC- (13) 

Notice the slight abuse of notation: is not the aa subma- 
trix of E. In fact El and Ej? are both matrices in the central 
region indices. Combining Eqs. (fTTI ) and ( fT2b we can write 

Gc{t,t') = 5o,c(t,t') 

+ dTidT2go,c{T, Ti)Etot(Ti, T2)Gc(r2, i<l)4) 



which expresses Gc in terms of the equilibrium propagator of 
the non-interacting, uncoupled system, go, and the total self- 
energy 

Et„t = E + El + E,j. (15) 

C. Real-time Green's functions 

In order to evaluate expectation values of single-particle 
observables we need the real-time correlation functions. We 
work with two correlation functions also called the lesser and 
greater GFs and defined as 

G<{t,t') - iTr{gmC^H,jit')cHAt)} (16) 
G>(i,<') = -zTr{ft„CH^,(t)4/t')}. (17) 

Two other important real-time GFs are the retarded and ad- 
vanced GFs defined by 

Gr,(t,0 - 0it-t'){G>{t,t')-G<it.t')) (18) 
G?,(t,0 = e{t' -t){G<{t,t')-G>it,t')). (19) 

The four GFs are related via 

G> - G< = G*^ - G". (20) 

The lesser and greater GFs are just special cases of the 
contour-ordered GF. For example G^ (t, t') — G{t, t') when 
T = t is on the upper branch of C and r' = t' is on the lower 
branch. This can be used to derive a set of rules, sometimes 
referred to as the Langreth rules, for converting expressions 
involving contour-ordered quantities into equivalent expres- 
sions involving real-time quantities. We shall not list the con- 
version rules here, but refer to Ref. 39 (no initial correlations) 
or Ref. [3^ (including initial correlations). The usual procedure 
in non-equilibrium is then to derive the relevant equations on 
the contour using the standard diagrammatic techniques, and 
subsequently converting these equations to real time by means 
of the Langreth rules. An example of this procedure is given 
in Sec. nil Bl where the non-equilibrium GW equations are de- 
rived. 



1. Equilibrium 

In equilibrium, the real-time GFs depend only on the time 
difference t' — t. Fourier transforming with respect to this 
time difference then brings out the spectral properties of the 
system. In particular the spectral function 

A{ijj) = i[G''(w) - G°H] = i[G>{ijj) ~ G<H] (21) 

shows peaks at the quasiparticle (QP) energies of the sys- 
tem. In equilibrium we furthermore have the fluctuation- 
dissipation theorem, 

G<(w) = if{uj - ii)A{u) (22) 
G>H = f{u^y))A{u), (23) 
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relating the correlation functions to the spectral function and 
the Fermi-Dirac distribution function, /. The fluctuation- 
dissipation theorem follows from the Lehman representation 
which no longer holds out of equilibrium, and as a conse- 
quence one has to work explicitly with the correlation func- 
tions in non-equilibrium situations. 

2. Non-equilibrium steady state 

We shall work under the assumption that in steady state, all 
the real-time GFs depend only on the time-difference t' — t. 
Taking the limit —oo this will allow us to use the Fourier 
transform to turn convolutions in real time into products in fre- 
quency space. Applying the Langreth conversion rules to the 
Dyson equation ( fT4l l. and Fourier transforming with respect to 
t' — t then leads to the following expression for the retarded 
GF of the central region 

GU^) = 9lc{^)+glc{^WtoMGh{^)- (24) 
This equation can be inverted to yield the closed form 

G^(w) = [{u + iT^)Ic -he- - - 

(25) 

The equation for is obtained by replacing r by a and 77 by 
—J] or, alternatively, from G" = (G"")^. For the lesser corre- 
lation function the conversion rules lead to the expression 

G</> = G^E</>G^(u;) + A</> (26) 

where 

A</> = [Ic + G^ELJffo^^ [Ic + SLG^]. (27) 

The w-dependence has been suppressed for notational sim- 
plicity. Using that E^/f = (5S/^)-1-(G;;/")-i together with 
the equilibrium relations = —f{uj—lic)[9oc~9oc] ^^'^ 
9o,c = ~if(^ - A*c) - l)[.95,c - 9o,cl we find that 

A<(co) = 2tf]f{i,-fic)Gh{i^)Gl.{uj) (28) 
A>H = 2i,7[/(co - ^c) - 1]G^HG^H. (29) 

If the product G''(ti')G''(a;) is independent of ry we can con- 
clude that A(aj) in the relevant limit of small 77. How- 
ever, as explained below, this is not always the case. 

3. Bound states and the A-term 

We first focus on non-interacting electrons. In this case the 
non-equilibrium correlation functions 5 ^ ^ ^ must be evaluated 
from Eq. (|26] l with Efot = '^l + '^r- For energies outside 
the band-width of the leads we have — EJ^ = such that 
no broadening of the (non-interacting) levels is introduced by 
the coupling to the leads. At such energies we have — 
Qq = 2irig^g^, and we conclude from Eqs. (|28]l,(|29Tl that 
A</> becomes proportional to the spectral function, A — 



9c ~ 9c- Since A{uj) does not necessarily vanish outside the 
band-width of the leads (it has delta peaks at the position of 
bound states), it follows that A</> should be included in the 
calculation of g^^^ to properly account for the bound states. 
It is interesting to notice that /ic which defines the initial state 
of the central region, drops out of the equations for g if and 
only if there are no bound states. 

When interactions are present in the central region corre- 
lation effects will reduce the lifetime of any single-particle 
state in C. Mathematically, this is expressed by the fact that 
j^r _ Y,a ^jjj non-zero for all physically relevant energies. 
Consequently, the product G^ {u!)G°'{uj) will approach a finite 
value as ?7 ^ leading to a vanishing A</>. 

In conclusion, the A terms of Eqs. (|28]l,(|29]l always vanish 
when interactions are present in G, while for non-interacting 
electrons they vanish everywhere except for w corresponding 
to bound states. We mention that it has recently been shown 
in the time-dependent NEGF framework that the presence of 
bound states can affect the long time behavior of the current 
in the non-interacting case^i. 

III. THE GW EQUATIONS 

In this section we derive and discuss the non-equilibrium 
GW and second order Born (2B) approximations. However, 
before addressing the expressions for the self-energies we in- 
troduce an effective interaction which leads to a particularly 
simple form of the equations and at the same time provides a 
means for reducing self-interaction errors in higher order dia- 
grammatic expansions. 

A. Effective interaction 

The direct use of the full interaction Eq. (O results in a four- 
index polarization a function. The numerical representation 
and storage of this frequency dependent four-index function is 
very demanding, and for this reason we consider the effective 
interaction defined by 

VeS= ^ Vi<y,ja'cl„c!j„,Cj„'Ci„, (30) 
ij, (T(t' 

where 

^icr.ja' — ^ij-ij ^aa'^ijji- (31) 

This expression follows by restricting the sum in the full in- 
teraction Eq. Q to terms of the form Vij^ijcl^Cj^,Cja'Cicr and 

The effective interaction is local in orbital space, i.e. it 
is a two-point function instead of a four-point function and 
thus resembles the real-space representation. Note, however, 
that in contrast to the real-space representation Via.ja' is spin- 
dependent. In particular the self-interactions, Via^ia, are zero 
by construction and consequently self-interaction (in the or- 
bital basis) is avoided to all orders in a perturbation expan- 
sion in powers of V. Since the off-diagonal elements (i ^ j) 
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of the exchange integrals Vij,ji are small, on expects that the 
main effect of the second term in Eq. dSTl ) is to cancel the 
self-interaction in the first term. 

It is not straightforward to anticipate the quality of a GW 
calculation based on the effective interaction ( l30l l as compared 
to the full interaction (O. Clearly, if we include all Feynman 
diagrams in E, we obtain the exact result when the full in- 
teraction (|3]l is used, while the use of the effective interaction 
dSOl l would yield an approximate result. The quality of this ap- 
proximate result would then depend on the basis set, becom- 
ing better the more localized the basis functions and equal to 
the exact result in the limit of completely localized delta func- 
tions where only the direct Coulomb integrals Vij^ij will be 
non-zero. 

However, when only a subset of all diagrams are included 
in E the situation is different: In the GW approximation only 
one diagram per order (in V) is included, and thus cancella- 
tion of self-interaction does not occur when the full interac- 
tion is used. On the other hand the effective interaction (ISTT i is 
self-interaction free (in the orbital basis) by construction. The 
situation can be understood by considering the lowest order 
case. There are only two first order diagrams - the Hartree 
and exchange diagrams - and each cancel the self-interaction 
in the other More generally, the presence of self-interaction in 
an incomplete perturbation expansion can be seen as a viola- 
tion of identities of the form (-Ic^g./ ■ ■ ■ Cia-Ci„ ■ ■ ■ Cja-"\-) = 0, 
when not all Wick contractions are evaluated. Such expecta- 
tion values will correctly vanish when the effective interaction 
is used because the prefactor of the CiaCic operator, Via.ia, is 
zero. The presence of self-interaction errors in (non-self con- 
sistent) GW calculations was recently studied for a hydrogen 



'^GW,i3{T-, t') 

W^{t,t') 

It is important to notice that in contrast to the conven- 
tional real-space formulation of the GW method, the spin- 
dependence cannot be neglected when the effective interac- 
tion is used. The reason for this is that V is spin-dependent 
and consequently the spin off-diagonal elements of W will 
influence the spin-diagonal elements of G, E, and P. A dia- 
grammatic representation of the GW approximation is shown 
in Fig. |2] 

As they stand, equations ([33]l-([35Tl involve quantities of the 
whole system (leads and central region). However, since Vij 
is non-zero only when i,j e C, it follows from Eq. ( l34l l. 
that W and hence E also have this structure. Consequently, 
the subscript C can be directly attached to each quantity in 
Eqs. ([33]l-(l35]l, however, for the sake of generality and nota- 
tional simplicity we shall not do so at this point. It is, however. 



atom . 

In App. |B] we compare the performance of the effective 
interaction with exact results for the Hartree and exchange 
self-energies of a benzene molecule. These first order results 
indicate that the accuracy of GW calculations based on the ef- 
fective interaction dSOl l should be comparable to GW calcula- 
tions based on the full interaction We stress, however, that 
in practice only the correlation part of the GW self-energy 
(second- and higher order terms) is evaluated using V^ff , while 
the Hartree and exchange self-energies are treated separately 
at a higher level of accuracy, see Sec. IV CI 



B. Non-equilibrium GW self-energy 

It is useful to split the full interaction self-energy into its 
Hartree and exchange-correlation parts 

E(r,r') = Eft(T,r') + E,,(T,r'). (32) 

The Hartree term is local in time and can be written 
E/j(r, r') — E/j(r)(5c(T, r') where Sq is a delta function 
on the Keldysh contour. Within the GW approximation 
the exchange-correlation term is written as a product of the 
Green's function, G, and the screened interaction, W, cal- 
culated in the random-phase approximation (RPA). With the 
effective interaction ( l30l l the screened interaction and the po- 
larization are reduced from four- to two-index functions. For 
notational simplicity we absorb the spin index into the orbital 
index, i.e. (ia) i (but we do not neglect it). The GW 
equations on the contour then read 

(33) 
(34) 

(35) 
I 

important to realize that the GF appearing in the GW equa- 
tions includes the self-energy due to the leads. 

Using the Langreth conversion rules^^ the retarded and 
lesser GW self-energies become (on the time axis), 

^GW,^J (t) = *G[, {t)W> (t) + iG< {t)w:^ it) (36) 
S^G^.W = *G</>(t)W,f >W, (37) 

where we have used the variable t instead of the time differ- 
ence t' — t. For the screened interaction we obtain (in fre- 
quency space), 

Wiuj) = V[I - P''iu;)V]-'^ (38) 
W</>{uj) = W{ij)P</> (1^)^(0;). (39) 



= iG^,{T,T'+)W,,{T,T') 

= Vrj 5c (r, r') + V / dn V^uPki (r, ti)Wi, (n , t') 
ki 

= -iGij{T,T')Gji{T' ,t). 
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where all quantities are matrices in the indices i, a and matrix 
multiplication is implied. Notice that the spin off-diagonal 
part of V will affect the spin-diagonal part of through the 
matrix inversion. 

Finally, the real-time components of the irreducible polar- 
ization become 

P^t) = -zG:^{t)G<{-t) - tG<{t)G%{-t) (40) 

i^f^W - ~^G</>{t)G>/<{-t). (41) 

From their definitions it is clear that both the polarization and 
the screened interaction obey the relations P/J- (w) — P^^{—lo) 
and W°j{u}) = Wj^{—u)), while for the self-energy and GFs 
we have E^^y(w) = '^^wi^Y and G"(w) = G"^(w)t. In 
addition all quantities fullfill the general identity X'' —X'^ = 
X^ — X"^. Notice that the GFs entering the GW equations are 
the We mention that equations similar to those derived above 
without the extra complication of coupling to external leads, 
have previously been used to calculate bulk bandstructures of 
excited GaAs>i^ 

In deriving Eqs. ( I38I39I I we have made use of the conver- 
sion rules dc'^i^^t') = and Sl^''{t,t') = 6{t - t'). With 
these definitions the applicability of Langeth rules can be ex- 
tended to functions containing delta functions on the contour. 
Notice, however, that with these definitions relation ( fTsT i does 
not hold for the delta function. The reason why the delta func- 
tion requires a separate treatment is that the Langreth rules are 
derived under the assumption that all functions on the contour 
are well behaved, e.g. not containing delta functions. 

We stress that no spin symmetry has been assumed in the 
above GW equations. Indeed by reintroducing the spin index, 
i.e. i — > [icr) and j — > (jcr'), it is clear that spin-polarized 
calculations can be performed by treating Gn and Gn inde- 
pendently. 

Within the GW approximation the full interaction self- 
energy is given by 

I](r,r') = S^(t,t') + W(r,r'), (42) 

where the GW self-energy can be further split into an ex- 
change and correlation part, 

W(t,t') = S,(5c(t,t') + S]con-(r,r'). (43) 

Due to the static nature of and Y^x we have 

^ !]</> ^ 0. (44) 

The retarded components of the Hartree and exchange self- 
energies become constant in frequency space, and we have 
(note that for Yih and Y.x we do not use the effective interac- 
tion ^) 

^M. = -^T.'^uii^mks (45) 

M 

= ^Y.Gti{i^mKi3- (46) 

kl 

Due to ( |44] |. it is clear that Eq. (|37| i yields the lesser/greater 
components of Scon-. Since Yj^ouij, r') does not contain delta 



functions its retarded component can be obtained from the re- 
lation, 

- Q{-t)\Y.l^{t) - Y.<y,(t)\. (47) 

The separate calculation of and S^q,,. from Eqs. (I46ll.(l47ll 
as opposed to calculating their sum directly from Eq. ( |36] |. 
has two advantages: (i) It allows us to treat E^r, which is the 
dominant contribution to Sgvi^, at a higher level of accuracy 
than Scon^ see App. |A] (ii) We avoid numerical operations 
involving C and W"^ in the time domain, see App.|E] 




FIG. 2: The GW and second Born self-energies, Egw and E2S, 
can be obtained as functional derivatives of their respective $- 
functionals, ^gwIG] and $2b[G]. Straight lines represent the full 
Green's function, G, i.e. the Green's function in the presence of 
coupling to the leads and interactions. Wiggly lines represent the 
interactions. 



C. Non-equilibrium second Born approximation 

When screening and/or strong correlation effects are less 
important, as e.g. in the case of small molecules, the higher- 
order terms of the GW approximation are small and it is more 
important to include all second order diagrams^^. The full sec- 
ond order approximation, often refeiTed to as the second Born 
approximation (2B), is shown diagrammatically in Fig.|2l As 
we will use the 2B for comparison with the GW results we 
state the relevant expressions here for completeness. The non- 
equilibrium 2B has recently been applied to study atoms in 
laser fields'*'*. 

On the contour the 2B self-energy reads (with the effective 
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interaction (|30] |) 

S2B,^,(r,r') = ^G,,(T,T')Gfci(T,r')G/fc(T',T)T/,feV;,, 

kl 

- G^k{T, T')Gki (r', t)Gi, (t, rOV^iV^jfe 

kl 

(48) 

Notice that the first term in is simply the second order 
term of the GW self-energy. From Eq. ( l48T l it is easy to obtain 
the lesser/greater self-energies, 

S2<£.w = T.Gf/''it)Gt/''{t)G>/<i~mkV,i 

kl 
kl 

where t has been used instead of the time difference t — t' . 
Since these second order contributions do not contain delta 
functions of the time variable, we can obtain the retarded self- 
energy directly from the Kramers-Kronig relation 

^2B{t) = e{^t)[j:>s{t) - i]<B(t)], (49) 

see App.|El 

IV. CURRENT FORMULA AND CHARGE 
CONSERVATION 

In this section we address the question of charge conserva- 
tion in the model introduced in Sec. Ill Al In particular, we 
ask under which conditions the current calculated at the left 
and right sides of the central region are equal, and we show in 
Sec. IIV E)l that this is fulfilled whenever the self-energy used 
to describe the interactions is <I>-derivable, independently of 
the applied basis set. 

A. Current formula 

As shown by Meir and Wingreen^, the particle current 
from lead a into the central region can be expressed as 

Ic^J ^Tr[l]<(c^)G>(c.) - J:>{cu)G<{co)], (50) 

where matrix multiplication is understood. By writing / = 
{II — Ir)/"^ one obtains a current expression symmetric in 
the L, R indices, 

^^ij MrL-rfl)G<+(/Lri-/flr^)(G^-G^)]dc. 

(51) 

where we have suppressed the uj dependence and introduced 
the coupling strength of lead a. Fa — — S°]. We note 
in passing that for non-interacting electrons the integral has 
weight only inside the bias window whereas this is no longer 
true when interactions are present. 



B. Charge conservation 

Due to charge conservation we expect that in steady-state 
~ —In — I, i.e. the current flowing from the left lead to 
the molecule is the negative of the current flowing from the 
right lead to the molecule. Below we derive a condition for 
this specific form of particle conservation. 

From Eq. dSOl l the difference between the currents at the left 
and right interface. A/ = /l + /_r, is given by 

A/ = 1 ^Tr[(S< + I]<)G> ~ (E> + S>)G<] (52) 

To obtain a condition for A/ = in terms of S we start by 
proving the general identity 

j ^Tr[E<,(c.)G>(c.) - ^>,{uj)G<{^)\ = 0. (53) 

To prove this, we insert G</> = Gj;i;J;{>Gg, + A</> (from 
Eq. ( |26l )) in the left hand side of Eq. ( |53] l. This results in two 

terms involving G^TifJf^ G"" and two terms involving A</>. 
The first two terms contribute by 

j ^Tr[l]<,G'^S>,G'^ - E>,G'^E<,G'^] . (54) 

Inserting E>4 = I]<t + (G")'i - {G'')~^ (see note SI) in 
this expression and using the cyclic invariance of the trace, it is 
straightforward to show that Eq. ( |54l i vanishes. The two terms 
involving A^/^ contribute to the left hand side of Eq. ( |53] | by 

j ^Tr[S<,(c.)A>(c.) - I]>,(c.)A<(a;)]. (55) 

As discussed in Sec. dll C 3l l A< and A> are always zero 
when interactions are present. In the case of non-interacting 
electrons we have = ^l^^ + ^i?^^' which vanish out- 

side the band width the leads. On the other hand A</> is 
only non-zero at energies corresponding to bound states, i.e. 
states lying outside the bands, and thus we conclude that the 
term ( |55] l is always zero. 

From Eqs. (|52] | and ( |53] | it then follows that 

A/ = j ^Tr[S<HG>H-I]>HG<H]. (56) 

We notice that without any interactions particle conservation 
in the sense A/ = is trivially fulfilled since S = 0. When 
interactions are present, particle conservation depends on the 
specific approximation used for the interaction self-energy, E. 

C. Conserving approximations 

A self-energy is called conserving, or ^-derivable, if it 
can be written as a functional derivative of a so-called $- 
functional, E[G] = 5^G]/5G?-^ Since a ^-derivable self- 
energy depends on G, the Dyson equation must be solved self- 
consistently. The resulting Green's function automatically ful- 
fills all important conservation laws including the continuity 
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equation which is of major relevance the context of quantum 
transport. 

The exact ^[G] can be obtained by summing over all skele- 
ton diagrams, i.e. closed diagrams with no self-energy inser- 
tions, constructed using the full G as propagator Practical ap- 
proximations are then obtained by including only a subset of 
skeleton diagrams. Two examples of such approximations are 
provided by the GW and second Born ^-functional and as- 
sociated self-energies which are illustrated in Fig.|2l Solving 
the Dyson equation self-consistently with one of these self- 
energies thus defines a conserving approximation in the sense 
of Baym. 

The validity of the conservation laws for <i>-derivable self- 
energies follows from the invariance of $ under certain trans- 
formations of the Green's function. For example it follows 
from the closed diagramatic structure of <1> that the transfor- 
mation^^ 

G(rT,r'T') e'^'^''^^GirT,r'T')e-'^'-'''^'\ (57) 

where A is any scalar function, leaves ^[G] unchanged. Us- 
ing the compact notation (ri,ri) = 1, the change in $ 
when the GF is changed by 6G can be written as 5$ = 
/ dld2I](l, 2)SG{2, 1+) = 0, where we have used that S = 
6^[G]/SG. To first order in A we then have 

S<^> ^ i J dld2I](l,2)[A(2) - A(1)]G'(2,1+) 

^ i J dld2[E(l,2)G(2,l+) -G(1,2+)S(2,1)]A(1). 

Since this hold for all A (by a scaling argument) we conclude 
that 

J d2[S(l, 2)G(2, 1+) - G(l, 2+)E(2, 1)] = 0. (58) 

It can be shown that this condition ensures the validity of the 
continuity equation (on the contour) at any point in space^i. 



which in matrix notation takes the form 
j dr'Tr[I](r,r')G(T',T+)-G(r-,r')S(T',T)] = 0. (60) 

Here is exactly the self-energy matrix obtained when the 
diagrams are evaluated using Gy and the Vij,ki from Eq. ((Sj. 
The left hand side of Eq. (|60l l, which is always zero for a 
derivable S, can be written as Tr[A< {t, t)] when A is given 
by Eq. dCTT i with B = E and G = G. It then follows from 
the general result (IC2l i and the condition ( |56] | that current con- 
servation in the sense II = —Ir is always obeyed when E is 
^-derivable. 

The above derivation of Eq. (|60] | relied on all the Coulomb 
matrix elements, Vijki, being included in the evaluation of E. 
Thus the proof does not carry through if a general truncation 
scheme for the interaction matrix is used. However, in the 
special case of a truncated interaction of the form dSOl l. i.e. 
when the interaction is a two-point function, Eq. ( l60l l remains 
valid. To show this, it is more appropriate to work entirely in 
the matrix representation and thus define <I>[G,;j (r, r')] as the 
sum of a set of skeleton diagrams evaluated directly in terms 
Gij and Vij . With the same argument as used in Eq. (|57] |. it 
follows that $ is invariant under the transformation 

G,,(r,r') ^ e*^'(^)G,j(T,T')e-^^^(^'\ (61) 

where A is now a discrete vector. By adapting the arguments 
following Eq. ( |57] i to the discrete case we arrive at Eq. ( |58] l 
with the replacements ri ^ i and Y2 j and with the inte- 
gral replaced by a discrete sum over j. Summing also over i 
leads directly to Eq. ( I6OI 1 which is the desired result. 

To summarize, we have shown that particle conservation 
in the sense II = —Ir, is obeyed whenever a ^-derivable 
self-energy is used and either (i) all Coulomb matrix elements 
VijM or (ii) the truncated two-point interaction of Eq. ( |30] |. 
are used to evaluate E. 



V. IMPLEMENTATION 



D. Charge conservation from ^-derivable self-energies 

Below we show that A/ of Eq. ( |56] l always vanishes when 
the self-energy is <I>-derivable, i.e. the general concept of a 
conserving approximation carries over to the discrete frame- 
work of our transport model. 

We start by noting that Eq. (l58T l holds for any pair 
G(l, 2), E[G(1, 2)] provided E is of the <I>-derivable form. In 
particular Eq. ( |58] l does not assume that the pair G, E[G] ful- 
fill a Dyson equation. Therefore, by taking any orthonormal, 
but not necessarily complete set, {0^}, and writing G(l, 2) — 
4'i{i^i)Gij{Ti,T2)(f>*{r2) we get from Eq. (|58ll after inte- 
grating over ri, 

^ / dr'[E,,(r,r')G,,(r',T+)-G.,(r^,r')E,,(r',T)]=0, 



In this section we describe the practical implementation of 
the Wannier-GVF transport scheme. After a brief sketch of the 
basic idea of the method we outline the calculation of the non- 
interacting Hamiltonian matrix elements and Coulomb inte- 
grals in terms of Wannier orbitals. The explicit expression for 
the Green's function is given in Sec. lVDl and in Sec. lVFI we 
describe our implementation of the Pulay mixing scheme for 
performing self-consistent Green's function calculations. We 
end the section with a discussion of the present limitations and 
future improvements of the method. 

A. Interactions in tlie central region 

Most first-principles calculations addressing transport in 
molecular contacts are based on the assumption that the 
charge carriers (electrons) can be considered as independent 
particles governed by an effective single-particle Hamiltonian. 
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A popular choice for the effective Hamiltonian is the Kohn- 
Sham Hamiltonian of DFT, 

hs ^ -^V^ + ue^t(r) + Vh(v) + v^^{v), (62) 

where Uext(r) is the external potential from the ions, u/i(r) 
is the classical Hartree field, and v-j-dv) is the exchange- 
correlation (xc-) potential which to some degree includes e-e 
interaction effect beyond the Hartree level. 

In the present method we rely on the KS Hamiltonian to de- 
scribe the metallic electrodes as well as the coupling into the 
central region, but replace the local xc -potential by a many- 
body self-energy inside the central region where correlation 
effects are expected to be most important. Clearly, this divi- 
sion does not treat all parts of the system on the same foot- 
ing, and one might be concerned that electrons can scatter off 
the artificial interface defined by the transition region between 
the mean-field and many-body description and thus introduce 
an artificial "contact resistance". Such unphysical scattering 
is certainly expected to affect the calculated properties if the 
transition region is very close to the constriction of the con- 
tact. On the other hand, the central region can, at least in 
principle, be chosen so large that the transition region occurs 
deep in the electrodes far away from the constriction. In this 
case the large number of available conductance channels in the 
electrodes should ensure that the calculated properties are not 
dominated by interface effects and the non-interacting part of 
the electrodes will mainly serve as particle reservoirs whose 
precise structure is unimportant. Thus the assumption of in- 
teractions in the central region seems justified in principle al- 
though it might be difficult to fully avoid artificial backscat- 
tering in practice. 

B. Wannier Hamiltonian and Coulomb integrals 

In order to make the evaluation and storing of the GW self- 
energy feasible, we use a minimal basis set consisting of max- 
imally localized, partially occupied Wannier functions'*^ ob- 
tained from the plane-wave pseudopotential code Dacapc'*^. 
Below we outline how the Hamiltonian is evaluated in the WF 
basis and refer to Ref. 49 for more details. 

The WFs used to describe the leads are obtained from a 
bulk calculation (or supercell calculation if the leads have fi- 
nite cross section). We define the extended central region (C2) 
as the molecule itself plus a portion of the leads. C2 should be 
so large that it comprises all perturbations in the KS potential 
arising from the presence of the molecular contact such that a 
smooth transition from C2 into the bulk is ensured. The WFs 
inside C2 are obtained from a DFT calculation with periodic 
boundary conditions imposed on the supercell containing C2. 
The resulting WFs will inherit the periodicity of the eigen- 
states, however, due to their localized nature they can be un- 
amigously extended into the lead regions. Thanks to the large 
size of C2, hybridization effects between the molecule and the 
metal leads will automatically be incorporated into the WFs. 
With the combined set of WFs (leadH-C2), we can then repre- 
sent any KS state of the contacted system up to a few electron 



volts above the Fermi energy. . 

In practice, the requirement of complete screening means 
that 3-4 atomic layers of the lead material must be included in 
C2 on both sides on the molecule. While this size of systems 
can be easily handled within DFT it may well exceed what is 
computationally feasible for a many-body treatment such as 
the GW method even with the minimal WF basis. For this 
reason we shall allow the central region (C) to consist of a 
proper subset of the WFs in C2, subject to the requirement 
that there is no direct coupling across it, i.e. (</),; | /is — 
for i ^ L and j E R where the left (right) lead by definition 
is all WFs to the left (right) of C. With this definition of C, 
the KS potential outside C is not necessarily periodic (this is, 
however, always the case outside C2), and consequently the 
calculation of the coupling self-energies becomes somewhat 
more involved as compared to the usual situation of periodic 
leads, see discussion in App. [D] We stress that the transmis- 
sion function for the non-interacting KS problem is exactly 
the same whether C or C2 is used as the central region as 
long as there is no direct coupling across region C. 

Having constructed the WFs we calculate the matrix ele- 
ments of the effective KS Hamiltonian of the contacted, unbi- 
ased system, {(l)i\hs\(t)j). To correct for double counting when 
the GW self-energy is added, we also need the matrix ele- 
ments, {4ii\vxc\<t>j), for WFs belonging to the central region. 

The matrix elements defining the interaction V in Eq. (O 
are calculated as the (unscreened) Coulomb integrals 

V.jM = J J drdr^ ^'^ > ' ^ ^ (63) 

for WFs belonging to the central region. The Coulomb inte- 
grals are evaluated in Fourier space using neutralizing Gaus- 
sian charge distributions to avoid contributions from the peri- 
odic images, see notei5(X 

C. Hartree and exchange 

As already mentioned it is not feasible to include all the 
interaction matrix elements when evaluating the frequency- 
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FIG. 3: The extended central region (C2) is chosen so large that 
it comprises all perturbations in the effective DFT potential arising 
from the molecular contact. The central region (C) can be a proper 
subregion of C2, but it must be so large that there is no direct cou- 
pling across it. We solve for the self-consistent Kohn-Sham potential 
within C2, but replace the static xc potential by the GW self-energy 
inside C. 
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dependent part of the many-body self-energy, Econ , which is 
therefore calculated using effective interaction of Eq. (|30l ). 

However, the exchange term, which can be unambigu- 
ously separated from the GW self-energy, is evaluated 
from Eq. ( |46] | using all Coulomb elements of the forms 
{Vijji}, {Viijj}, {Viijj}}. As shown in appendix 
|A]this produces results within 5% of the exact values. 

The KS Hamiltonian already includes the Hartree potential 
of the DFT groundstate. In a self-consistent, finite bias GW 
calculation the relevant Hartree potential will deviate from the 
DFT Hartree potential due to the finite bias and the fact that 
the xc -potential is replaced by the GW self-energy. This cor- 
rection, which is much smaller than the full Hartree potential. 



Several comments are in order. First, we notice that all quanti- 
ties except for v^c, hg, and E)^ [Ssf*'], are bias-dependent, how- 
ever, to keep the notation as simple as possible we omit any 
reference to this dependence. The terms E£ and E^ account 
for the coupling to the leads. By subtracting v^c from hg 
we ensure that exchange-correlation effects are not counted 
twice when we add the GW self-energy, '^^q-^- The term 
Awft = 5][JG] — EJ^ [gi''''*] is the change in Hartree poten- 
tial relative to the equilibrium DFT value. This change is due 
to the applied bias and the replacement of Vxc by E^^y (even 
in equilibrium the Hartree field will change during the GW 
self-consistency cycle). The Hartree potential in C originat- 
ing from the electron density in the electrodes, which enters 
G^ through hg, is assumed to stay constant when the system 
is driven out of equilibrium, i.e. the out-of-equilibrium charge 
distribution in the leads is assumed to equal the equilibrium 
one. 

Finally, in order to make contact with the general formalism 
of Sec. (HJl, and in particular Eq. (IZST i, we note that the ma- 
trix elements hij defining the effective single-particle Hamil- 
tonian in Eq. ([T]), are related to the quantities introduced above 
via 

{{(t)i\hg - Vxc\4>]) - ^hld's'^Vj fori, j both in C 
{(j)i\hg\4)j) + {^lL{R) ~£F)Sij fori, j both in 
{4>i\hg\(j)j) otherwise 

E. Frequency dependence 

To represent the temporal dependence of the Green's func- 
tions and GW self-energies we use an equidistant frequency 
grid with Ng grid points and grid spacing 5. Thus the GFs 
(and the GW self-energies) are represented by TV^y x x iV^ 
matrices. At each of the discrete frequencies u!i = riid. 
Hi = . . . Ng, we have an iV^, x iV^, matrix representation 



is treated in the same way as the exchange term, i.e. calcu- 
lated from Eq. (l46T l with all Coulomb elements of the form 

{{VijAj},{Vi3,ji}A'^ii.3j}A^ii,ij}}- ^s for the exchange 
terms this yields results within 5% of the exact values, see 

El 



D. Expression for C 

To simplify the notation in the following we omit the sub- 
script C as all quantities will be matrices in the central region. 
The retarded GF of the central region is obtained from 



- (E[, [G] - E[, [.g?^>] ) - i:^aw [G]]-' (64) 
I 

of G{uJi) in the WF basis. The grid spacing, S, should be 
small enough that all features in the frequency dependence of 
the GFs and self-energies can be resolved. At the same time 
the frequency grid should be large enough (contain enough 
points) to properly describe asymptotic behavior (the tail) of 
the GFs. Although the tail is irrelevant for the current in 
Eq. ( ISTI l. it contributes to the self-energy, I^gw [G] . In prac- 
tice, Ng and S should be increased, respectively decreased, 
until the results do not change. 

To avoid time consuming convolutions on the frequency 
grid, we use the Fast Fourier Transform (FFT) to switch be- 
tween frequency and time domains. An important but tech- 
nical issue concerning the evaluation of retarded functions is 
discussed in App.|E] 

F. Self-consistency 

Since E depends on G, and G depends on E, the Dyson 
equations Eqs. ( |26l ) and ( |64l ) must be solved self-consistently 
in conjunction with the equations for the GW, Hartree, and 
exchange self-energies. In practice this self-consistent prob- 
lem is solved by iteration. Clearly, the iterative approach re- 
lies on the assumption that the problem has a unique solution 
and that the iterative process converges to this solution. For 
all applications we have studied so far this has been the case. 
In order to stabilize the iterative procedure, we use the Pu- 
lay scheme^ ' to mix the GFs of the previous N iterations very 
analogue to what is done for the electron density in many DFT 
codes. More specifically the input GF at iteration n is obtained 
according to 

j—n — N j—n — N 

(65) 

To determine the optimal values for the expansion coeffi- 
cients, c", we first define an inner product in the space of 
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(retarded) GFs 

(G-.^ G"-'^) f in,[GZiu^)rMG:;iiu;)]du;. (66) 

Equivalent inner products can be obtained e.g. by using the 
real part of the GF instead of the imaginary part or the lesser 
component instead of the retarded. The Pulay residue matrix 
determining the coefficients c" is then given by 

— ~ G^ouii Gl^^ — Goy-'t), (67) 

where i,j = n — N , . . . ~ 1. We typically use a mix- 
ing factor around a k, 0.4. During the mixing procedure one 
must keep track of both the retarded and lesser GF since one 
does not follow directly from the other However, it is im- 
portant that the same coefficients, c", are used for mixing the 
two components. If separate coefficients are used for and 
, the fundamental relation (|20] | is not guaranteed during the 
self-consistent cycle. As noted above we define the residue 
exclusively from the retarded GF. In practice we always find 
that once the retarded GF has converged, the lesser GF has 
converged too, and this justifies the use of common expansion 
coefficients for the two GF components. 



use Hartree-Fock or some other mean-field orbitals. Out of 
equilibrium the WFs will be distorted due to the change in 
electrostatic potential, however, this effect is not included. Al- 
though the manifold spanned by the WFs, i.e. the KS eigen- 
states up to a few electron volts above the Fermi level, are 
expected to represent the GW quasiparticle wave functions of 
the same energy range quite well, an accurate representation 
of the screened interaction might require inclusion of high- 
energy eigenstates. 

With the present implementation of the GW scheme it is 
not feasible to include more than a few electrode atoms in 
addition to the molecule itself in the GW region (region G in 
Fig.O. The use of a small central region region might affect 
the description of image charge formations in the electrode, 
and it might introduce artificial backscattering at the DFT-GW 
interface. 

The use of larger and more accurate basis sets as well the 
inclusion of more electrode atoms in the GW region are not 
fundamental but practical limitations of the method, which in 
principle could be removed by invoking efficient simplifica- 
tions/approximations into the present formalism. 



VI. ANDERSON MODEL 



G. Overview 

Below we give an overview of the various steps involved in 
performing a self-consistent non-equilibrium GW transport 
calculation: 

• Perform DFT calculations for the electrodes and the ex- 
tended central region (region C2 in Fig.|3]l. 

• Construct the Wannier functions, and obtain the matrix 
representation of the KS Hamiltonian for the contacted 
system in equilibrium. Evaluate the matrix elements for 
Vxc and relevant Coulomb integrals for Wannier func- 
tions belonging to the central region (C). 

• Fix the bias voltage, and calculate the coupling self- 
energies Eq. ( fTsl l as described in App. |D] (these stay 
unchanged during self-consistency). 

• Evaluate the initial (non-interacting) Green's functions, 

and G^, e.g. from the KS Hamiltonian. 

• From G^ and G^, construct the desired interaction 

self-energies (E/j, Y.^, Y,gw, or S2b)- 

• Test for self-consistency. In the negative, obtain a new 
set of output Green's functions from Eqs. ( l64l l and (|26] | 
and mix with the previous GFs as described in Sec. IV Fl 

H. Limitations and future improvements 

The main approximation of the present implementation is 
the use of a fixed, minimal basis set. We have used WFs ob- 
tained from the DFT-PBE orbitals, however, one could also 



Since its introduction in 1961 the Anderson impurity 
modeP^ has become a standard tool to investigate strong cor- 
relation phenomena such as local moments formation, Kondo 
effects and Coulomb blockade. The Anderson model de- 
scribes a localized electronic level of energy £c and correlation 
energy U coupled to a continuum of states. Thus the central 
region-part of the Hamiltonian reads 

He = £cC^c + Un^ni. (68) 

In equilibrium, accurate results for the thermodynamic 
properties of the Anderson model have been obtained from 
the Bethe ansatz^^''^'*, quantum monte carlo simulation a^^'^^ , 
and numerical renormalization group theory 

Out of equilibrium, the low-temperature properties of the 
Anderson model have been much less studied. The earliest 
work addressed the problem by applying second-order pertur- 
bation theory in the interaction strength U Despite the 
simplicity of this approach it provides a surprisingly good de- 
scription of the (equilibrium) spectral function. There are, 
however, several fundamental problems related to the non-self 
consistent low-order perturbative approach: (i) the result de- 
pends on the starting point around which the perturbation is 
applied, (ii) it inevitably violates the conservation laws, and 
(iii) it applies only in the small- [/ limit. Methods relying on 
the slave-boson technique^" have been developed to explore 
the strong correlation regime of the model. The noncrossing 
approximation is believed to work well in the infinite-C/ limit 
and for sufficiently small tunneling strength, F, but it fails to 
reproduce the correct Fermi liquid behavior at low temper- 
atures. '''■''^ More recently, a finite-f/ slave-boson mean-field 
approach^^ has been proposed. Finally we mention that a 
number of more advanced schemes have been used to address 
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non-equilibrium Kondo-like phenomea focusing on the low- 
energy properties of the Anderson model in the limit where U 
is much larger than the hybridization energy, F^i^^i^^ . 

While the Anderson model is normally used to describe 
strongly correlated systems, the main application of the GW 
approximation has been to weakly interacting quasi-particles 
in closed shell systems such as molecules, insulators and 
semi-conductors. In view of this, one could argue that the 
GW method is inappropriate for the Anderson model. Nev- 
ertheless, we find this application rather instructive as it illus- 
trates some general features of the GW approximation includ- 
ing the role of self-consistency both in relation to charge con- 
servation and the line shape of spectral functions. Moreover, 
as many important transport phenomena, like Kondo effects 
and Coulomb blockade, are well described by the Anderson 
model, it should always be of interest to benchmark a trans- 
port scheme against this model. 

In a very recent study the GW approximation was ap- 
plied to the Anderson model in equilibrium for interaction 
strengths U/T up to 8.4/0.65 « 13, and various tempera- 
tures. For the largest interaction strength it was found that 
GW prefers to break the spin symmetry leading to directly 
erroneous results in the Kondo regime. For intermediate in- 
teraction strengths (U/T — 4.2/0.65 w 6.5) where GW does 
not break the spin symmetry, it was concluded that GW does 
not describe the T-dependence of the Kondo effect well. Nev- 
ertheless, we show here that at T = the width of the GW 
Kondo-like resonance follow the analytical result for Tk quite 
well for intermediate interaction strengths. 

Here, as in our previous paper^^, we focus on the zero tem- 
perature, non-equilibrium situation. We consider interaction 
strengths of U/T up to 8 (we keep fix U — 4 and vary F). 
For these interaction strengths we always find a stable non- 
magnetic GW solution, i.e. G|| = Gn- In contrast, the HF 
solution can develop a magnetic moment for U/T > tt (de- 
pending on bias voltage and Ec). We adopt the wide-band ap- 
proximation where the coupling to the continuum is modeled 
by constant imaginary self-energies Sl + Sfl = —iT. With- 
out loss of generality we set Ep = 0. In all calculations the 
frequency grid extends from -15 to 15 with the grid spacing 
ranging from 0.1 to 0.0005. 



A. Equilibrium spectral function 

In Fig. |4] we show the ec-dependence of the equilibrium 
spectral function, A{uj) = — ImG"'(a;), for U — A and F = 
0.65. The HF solutions are Lorentzians centered at Ehf = 
£c + U{ncr) with a full width at half maximum (FWHM) given 
by 2F. As can be seen the position of the HF peaks do not 
vary linearly with Ec- Instead there is a "charging resistance" 
for the peak to move through the Fermi level due to the cost 
in Hartree energy associated with the filling of the level. This 
effectively pins the level to Ep. 

Moving from HF to the second Born approximation, the 
Lorentzian shape of the spectral peak is distorted due to the uj- 
dependence of the 2B self-energy. We can observe a general 
shift of spectral weight towards the chemical potential as well 



as a narrowing of the resonance as it comes closer to Ep. 

The redistribution of the spectral weight towards the chem- 
ical potential becomes even more pronounced in the GW ap- 
proximation. For F — [/ <ec<~r (the so-called Kondo 
regime) a sharp peak develops at Ep. For U/T sufficiently 
large the Kondo effect should reveal itself as a peak in the 
spectral function with a full width at half maximum (FWHM) 




FIG. 4: (color online). Spectral function of the central site for F — 
0.65, U = 4.0 and different values of £c. The inset in the lower 
panel is a zoom of the GW spectral peak around u) = 0. 
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FIG. 5: (color online). Full width at half maximum (FWHM) of the 
Kondo resonance as calculated in the GW approximation and from 
the analytical result Eq. ( |69b . The interaction strength isU — 4 and 
Ec is varied in the Kondo regime. 
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given approximately by the Kondo temperaturoSa 

Tk ^ 0.5(2rj7)i/2gxp[7rec(e^ + U)/2TU]. (69) 

In Fig. |5] we compare the above expression for Tk with the 
FWHM of the GW Kondo peak. The exponential scaling of 
Tk is surprisingly well reproduced. Deviations from the ex- 
ponential scaling naturally occur for smaller values of U /T 
(not shown) where the Kondo effect does not occur and ( |69] l 
does not apply. In accordance with recent work^, we were 
not able to obtain non-magnetic GW solutions in the strong 
interaction regime (U /T > 8). 

In Fig.|6]we show the dependence of the spectral function 
on the ratio U/T for the central level at the symmetric posi- 
tion Ec = —U /2 = —2. For U/r — 2 there is no significant 
difference between the three descriptions. This is to be ex- 
pected since the correlation plays a minor role compared to 
the hybridization effects. In the weakly coupled limit, how- 
ever, correlations become significant and as a consequence the 
2B and GW results changes markedly from the Lorentzian 
shape and show a Kondo-like peak at the metal Fermi level. 
The 2B approximation significantly overestimates the width 
of the Kondo peak, indicating, as expected, that the higher 
order RPA terms enhance the strong correlation features. 

For large U /T it is known^*'^- that the spectral function, 
in addition to the Kondo peak, should develop peaks at the 
atomic levels £c and Ec + U. We find that the self-consistent 
2B and GW approximations always fail to capture these side- 
bands and instead distribute the spectral weight as a broad 
slowly decaying tail. These findings agree well with previ- 
ous results obtained with the fluctuation-exchange approxi- 
mation^^, and with GW studies of the homogeneous electron 
gas which showed that self-consistency in the GW self-energy 
washed out the satelite structure in the spectrum^S. 



B. Non-equilibrium transport 

We now move to the non-equilibrium case and introduce 
a difference in the chemical potentials of the two leads. In 
Fig.|7]we show the zero-temperature differential conductance 
under a symmetric bias, /x^ //^ = ±V/2, as a function of Ec for 
U — 4 and F = 0.65. The dl/dV at bias voltage V has been 
calculated as a finite difference between the currents obtained 
from Eq. ( BTI ) for bias voltages V and V + 6V, respectively. 
The 2B result falls in between the HF and GW results, and 
for this reason we will focus on the latter two in the following 
discussion. 

For V — there is only little difference between the three 
results which all show a broad conductance peak reaching 
the unitary limit at the symmetric point Ec = —U/2. The 
physical origin of the conductance trace is, however, very dif- 
ferent: While the HF result is produced by coherent trans- 
port through a broad spectral peak moving rigidly through the 
Fermi level, the GW result is due to transport through a nar- 
row Kondo peak which is always on resonance (for Ec in the 
Kondo regime). In all cases the width of the dl/dV curve is 
approximately U. In the GW case this is because the Kondo 
peak develops only when the central level is half occupied, 
i.e. —U ^ Ec ^ 0. In HF, on the other hand, the dl/dV peak 
acquires a width on the order of U due to the charge pinning 
effect discussed in Sec. lVI Al 

The difference in the mechanisms leading to the HF and 
GW results is brought out clearly as V is increased: for 
V <^r the bias has little effect on the HF conductance while 
the GW conductance drops dramatically already at biases 
comparable to Tk due to suppression of the Kondo resonance 
at finite bias. The suppression of the Kondo resonance is due 
to quasi-particle (QP) scattering. While QP scattering does 
not affect the life-time of QPs at Ep in equilibrium, it does 
so at finite bias where ImEGwiEp) becomes non-zero. We 
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mention that we do not observe a splitting of the GW Kondo 
resonance at finite V^^. 

The peaks appearing in the d//dy at the largest bias {V = 
4) occur when the central level is aligned with either the lower 
or upper edge of the bias window. It is worth noticing that 
the height of these peaks are smaller than the value of IGo 
expected from on-resonant transport through a single level. 
The reason for this is two-fold: (i) The bias window only hits 
the resonance with one edge (either upper or lower edge), and 
consequently only half the spectral weight enters the bias win- 
dow when the voltage is increased by AV as compared to the 
low-bias situation, (ii) The self-consistent charging resistance 
discussed in Sec lVI Al pins the level to the edge of the bias 
window making the resonance follow the bias. 



C. The Go Wo approximation 

Non self-consistent, or one-shot, GW calculations can be 
performed by evaluating the screened interaction and GW 
self-energy from some trial non-interacting Green's function, 
Gq. The resulting GqWq approximation, with Go obtained 
from an LDA/GGA calculation, has been found to yield very 
satisfactory results for the band gaps of insulators and semi- 
conductors^'-^^. For this reason, and due to its significantly 
lower computational cost, this GqWq approach has generally 
been preferred over the self-consistent GW. One rather unsat- 
isfactory feature of the perturbative GoWo method is its Gq- 
dependence. However, as will be demonstrated below, a just 
as critical problem in non-equilibrium situations is its non- 
conserving nature. 

Before we apply the GaWa approximation to the Ander- 
son model, we need to address a certain issue which unfortu- 
nately has led to an error in our previous paper Ref. |27l (all 
conclusions from that paper are, however, unaffected by the 
mistake.) 



1. Instability of the non-magnetic ground state 

Consider a system which admits a spin polarized ground- 
state at the Hartree level (notice that Hartree and HF is equiv- 
alent for the Anderson model when the effective interaction of 
Eq. ( [3TI 1 is used), and let Go denote the GF obtained from spin 
unpolarized Hartree calculation. It turns out that the analyti- 
cal properties of the screenined interaction, [Go], evalu- 
ated from Go will be wrong. In particular W(5'[Go] will not be 
retarded as it should be. The reason is that the RPA response 
function is ill defined around the non-magnetic, and thus un- 
stable, Gq. The problem has been previously mentioned by 
J. A. Whita^, and was brought to the authors attention by C. 
Spataru. 

For certain parameter values, the HF groundstate of the 
Anderson model develops a finite magnetic moment. As a 
consequence the analytic properties of Wq as calculated from 
Eq. ( [38] l with the unpolarized Ghf become wrong. In our pre- 
vious paper Ref. 27, this problem was not recognized because 
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FIG. 8: (color online). Differential conductance as a function of ap- 
plied bias for (7 = 4, F = 0.65 and Sc = —4. For these parameters, 
the non-magnetic HF solution is stable for bias voltages smaller than 
^ 1.6. The GqWo approximation yields different currents at the 
left and right interfaces (AJ 7^ 0), and yields negative differential 
conductance at finite bias. 



we, for numerical efficiency, applied the Kramers-Kronig re- 
lation ( |47] | to obtain from I]< — S^, instead of using 
Eq. ( [36] l. Thus by construction our Y,^ was retarded. Specif- 
ically, this implies that the GqWq spectral function plotted in 
Fig. 1 of that paper, as well as the dl/dV curves in the middle 
panel of Fig. 2 for Ec in the interval —3.6 to —0.4, are in- 
correct. In fact there exists no non-magnetic GoT4^o[Ghf] so- 
lution in these cases. We stress, however, that all conclusions 
from our paper are unaffected by this mistake. In particular we 
show below that for parameter values leading to a stable non- 
magnetic HF groundstate, the Go Wo approximation still vio- 
lates charge conservation and gives unphysical results such as 
negative differential conductance. Moreover, we arrive at the 
same conclusions for Go Wo self-energies constructed from 
the spin polarized HF Green's function, in which case the in- 
stability problem does not occur at all. 



2. Results of the Go Wo approximation 

In Fig. [8] we show the calculated dl/dV for the Anderson 
model with T = 0.65 and Sc = -4 for the HF, GW, and 
Go Wo [Ghf] approximations. For these parameters, the non- 
magnetic HF solution is stable for bias voltages smaller than 
^ 1.6, such that the Go Wo approximation based on a non- 
magnetic Ghf is indeed meaningful in this parameter range. 
The Go Wo conductance has been obtained as a finite differ- 
ence between the currents obtained from Green's functions 
with self-energies Sgvi/'[Ghf(^)] and T,gw[Ghf{V + SV)], 
respectively, where Ghf(^) is the HF Green's function eval- 
uated self-consistently under a bias voltage V. 

From Fig. [8] we conclude that the GqWo approximation 
leads to unphysical results in the form of strong negative dif- 
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FIG. 9: (color online). Upper panel: Occupation of the central 
site as function of £c for U = A, T = 0.65, and bias V = 0.8. 
Notice that the HF solution breaks the spin symmetry for some Sc 
values. Middle panel: Current calculated in self-consistent GW 
and GoWoiGnF,!, Ghf,;]. Lower panel: Violation of the continu- 
ity equation measured as the difference between the currents in the 
left and right leads. 



ferential conductance. Moreover, as shown in the lower panel 
of the figure, the Go Wo approach gives different values for 
II and Ip. We note in passing that this symmetry break 
comes from the different chemical potentials of the left and 
right leads. Finally, we mention that the increasing behavior 
of AI/I as function of bias voltage seems to be a general ef- 
fect. 

As already mentioned the HF solution breaks the spin sym- 
metry for certain parameter values. Meaningfull GqWq re- 
sults can still be obtained in this case provided the self-energy 
is constructed from the spin polarized HF Green's function. 
Figs. |9l and [TOl compare the result of such calculations with 
self-consistent GW for two different values of the bias volt- 
age. From the figures we draw the conclusions: (i) The GqWo 
and GW currents agree when the level is alomst empty/filled 
(ii) The current calculated in GqWq show unphysical behavior 
in and close to magnetic regime (iii) The violation of charge 
conservation in GqWq is more severe when the current is 
large. 



VII. BENZENE JUNCTION 

In this section we apply the Wannier-GVF method to a more 
realistic nano junction, namely a benzene molecule coupled to 
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FIG. 10: (color online). Same as Fig.[9l but for bias voltage V = 4.0. 



featureless leads. In contrast to the Anderson model consid- 
ered in the preceding section, the benzene junction represents 
a closed-shell system with the Fermi level lying within the 
HOMO-LUMO gap leading to rather low transmission for all 
but the strongest molecule-lead coupling strengths. 

The use of featureless (wide-band) electrodes is convenient 
as it allows us to isolate the effects of the electron-electron 
interactions. The use of more realistic contacts with energy 
dependent spectral features would lead to an additional renor- 
malization of the molecular levels making a clear separation 
between xc- and contact effects more difficult. We stress, 
however, that the contacts only enter the theory through the 
coupling self-energies which can be calculated once and for 
all as in the standard NEGF-DFT approach. Thus the use of 
more realistic contact self-energies is straightforward. 

To describe the benzene molecule we first perform a DFT 
calculation for the isolated molecule, see note 74. The 
KS eigenstates are then transformed into maximally localized 
WFs, and the KS Hamiltonian and Coulomb integrals are eval- 
uated in the WF basis. For the interactions we use the trunca- 
tion scheme V'-^'> defined in App. |A]to evaluate Hartree and 
exchange self-energies. As shown in table U this leads to re- 
sults within ~ 5% of the exact values. We use the effective 
interaction Eq. ( [3T| l for the correlation part of the GW self- 
energy. In all calculations we have applied a frequency grid 
extending from -100 to 100 eV, and grid spacings in the range 
0.2 to 0.02, depending on the value of F. 

In Sec. I VII Al we show that the experimental ionization po- 
tential of the isolated benzene molecule is very well repro- 
duced with our GW scheme. In Sec. IVIIBI we investigate 
the role of the coupling strength, F, on the spectrum of the 
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benzene junction. Finally, in Sec. lVllCl we calculate the non- 
equilibrium conductance of the junction and compare various 
approximations for the xc self-energy. 



A. Spectrum of isolated benzene 

Within our general transport formalism we model the sit- 
uation of a free molecule by using a very weak coupling 
to the wide-band leads, see Fig fTTT a). The contacts merely 
act as particle reservoirs fixing the number of electrons on 
the molecule and providing an insignificant broadening (F = 
O.OSeV) of the discrete energy levels. We fix the Fermi lev- 
els of the electrodes to Ep = —3 eV which is approximately 
half-way between the HOMO and LUMO levels (the precise 
position of Ep within the gap is unimportant for the results 
presented in this section). 

In Fig.[l2]we show the total density of states (DOS), 



(70) 



, where the sum runs over all WFs on the molecule. We use 
three different approximations: (i) DFT-PBE (ii) Hartree-Fock 
(iii) fully self-consistent GW . We stress that our calculations 
include the full dynamical dependence of the GW self-energy 
as well as all off-diagonal elements. Thus no analytic ex- 
tension is performed, and we do not linearize the self-energy 
around the DFT eigenvalues to obtain an approximate quasi- 
particle equation as is done in standard GW calculations. 



a) 




b) 
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The spectral peaks seen in Fig. [12] occurring above (below) 
the Fermi level correspond to electron addition (removal) en- 
ergies. In particular, the HOMO level should coincide with the 
(vertical) ionization energy of the isolated molecule, which in 
the case of benzene is /gxp = —9.2 eV^'. The PBE functional 
overestimates this value by 3 eV, giving Jpbe — — 6.2 eV in 
good agreement with previous calculations-'". The HF and 
GW calculation yields /hf — —9.7 eV and /gvk = —9.3 eV, 
respectively. Because of the limited size of the WF basis. 
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FIG. 12: Density of states for a benzene molecule weakly coupled to 
featureless leads (F = 0.05). The common Fermi levels of the leads 
is indicated. Notice the characteristic opening of the band gap when 
going from DFT-PBE to HF, and the subsequent (slight) reduction 
when correlations are included at the GW level. 
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FIG. 11: (color online), (a) Illustration of a benzene molecule cou- 
pled to featureless electrodes with different chemical potentials, (b) 
Iso-surfaces for the 18 partially occupied Wannier functions used as 
basis functions in the calculations. The WFs are linear combina- 
tions of Kohn-Sham eigenstates obtained from a DFT-PBE plane- 
wave calculation. 



FIG. 13: The HF and GW HOMO-LUMO gap of the benzene 
molecule as a function of the coupling strength V. The difference 
between the curves represents the reduction in the gap due to the 
correlation part of the GW self-energy. This value increases with 
the coupling strength, as screening by electrons in the leads becomes 
more effective. 
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the perfect agreement between the GW and experimental val- 
ues should not be taken too strict. Indeed, more accurate HF 
calculations predicts a HOMO level around —9.2 eV which 
is 0.5 eV higher than our HF result. The deviation of our 
HF calculation from this number is two-fold: (i) The use of 
the truncated interaction V*^^) to evaluate the exchange self- 
energy introduces an error of ^ 0.1 eV, see table |I1 (ii) The 
difference between the PBE orbitals (from which our WFs are 
constructed) and the true HF orbitals. 

Returning to Fig. [121 we notice a dramatic opening of 
the HOMO-LUMO gap when going from PBE to HF (and 
GW). This effect is due to the inabiUty of the LDA/GGA 
functionals to fully cancel the spurious self-interaction con- 
tained in the Hartree potential. For the same reason, the self- 
interaction free HF method generally yields better spectra than 
the LDA/GGA functionals for small, localized systems where 
self-interaction terms are significant and dynamic screening is 
small. The GW spectrum resembles the HF spectrum with a 
slight reduction of the gap by ~ 1.0 eV. As we show in the 
next section, the GW gap shrinks as the coupling strength, F, 
is increased. 
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FIG. 14: (color online). Differential conductance of the benzene 
junction for Fl = Tr = 0.25 eV. Notice the strong Go dependence 
of the GoWo result. 



B. Contact enhanced screening (the role of F) 

In Fig. [13] we plot the size of the HOMO-LUMO gap as 
a function of the coupling strength F. Both the HF and GW 
gaps decrease as T is increased. For the HF gap, this is a 
simple consequence of the redistribution of charge from the 
HOMO to the LUMO when the resonances broaden and their 
tails start to cross the Fermi level. As this happens the HOMO 
(LUMO) self-interaction term in T^x will become less (more) 
negative, and consequently the HF gap shrinks. 

The GW quasi-particle energies consist of a HF eigenvalue 
and a correlation contribution coming from the real part of the 
dynamic GW self-energy. 

According to Fig.[T3j the correlation part of the QP gap. 
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increases significantly with F. In fact for a large range of cou- 
pling strengths, the reduction of the gap is more than 3 eV. 
This reduction can be understood from the enhanced mobility 
of the electrons on the molecule when the coupling is strong. 
The enhanced mobility allows for more efficient screening and 
this reduces the QP gap. The difference between the large- and 
small r limits is analogue to the difference between extended 
and confined systems. In extended systems where screening 
is significant, band gaps are overestimated by HF, and cor- 
relation contributions to the gap are large. In confined sys- 
tems, such as atoms and small molecules, screening effects 
are unimportant and HF usually yields good HOMO-LUMO 
gaps. 
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FIG. 15: (color online). Equilibrium DOS for the benzene molecule 
coupled to wide-band leads with a coupling strength of Tl ~ Fu = 
0.25 eV. Upper panel shows DFT-PBE and HF single-particle ap- 
proximations while the lower panel shows the self-consistent GW 
result as well as one-shot GqVKo results based on the DFT and HF 
Green's functions, respectively. 



C. Conductance 

In this section we consider the transport properties of the 
benzene junction under a symmetric bias, /iL/fl = ±V^/2, and 
a coupling strength of F^ = F/? = 0.25 eV. 

In Fig.[T4]we compare the differential conductance, d//dV^, 
calculated in self-consistent DFT-PBE, HF, GW, as well as 
non self-consistent GqWq using either the DFT-PBE or HF 
Green's function as Go. The d//dT^ has been obtained by 
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numerical differentiation of the I{V) curves calculated from 
Eq. ( ISTl ). For the DFT calculation the finite-bias effects have 
been included at the Hartree level, i.e. changes in the xc- 
potential have been neglected. We notice that the HP and 
Go Wo [Ghf] results are close to the self-consistent GW result. 
These approximations all yield a nearly linear IV with a con- 
ductance of ^ O.OSGo. In contrast the DFT and Go Wo [Gdft] 
yield significantly larger conductances which increase with 
the bias voltage. We note that the violation of charge conser- 
vation in the Go Wo calculations is not too large in the present 
case (A/// < 5%). This is in line with our general obser- 
vation, e.g. from the Anderson model, that AI/I grows with 
I. 

The trends in conductance can be understood by consider- 
ing the (equilibrium) DOS of the junction shown in Fig. [15] 
As for the free benzene molecule (see Fig. [T2] | the DFT 
HOMO-LUMO gap is much smaller than the HF gap, and 
this explains the lower conductance found in the latter case. 
The GW gap falls in between the DFT and HF gaps, how- 
ever, the magnitude of the DOS at Ep is very similar in GW 
and HF which is the reason for the similar conductances. It 
is interesting to notice that the HOMO-LUMO gap obtained 
in the Go Wo calculations resemble the gap obtained from 
Go, and that the self-consistent GW gap lies in between the 
Go Wo [Gdft] and GoWo[Ghf] gaps. 

The increase in the Go Wo [Gdft] conductance as a function 
of bias occurs because the LUMO of the Go Wo [Gdft] calcu- 
lation moves downwards into the bias window and becomes 
partly filled as the voltage is raised. In a self-consistent calcu- 
lation this would lead to an increase in Hartree potential which 
would in turn raise the energy of the level. The latter effect is 
missing in the perturbative Go Wo approach and this can lead 
to uncontrolled changes in the occupations as the present ex- 
ample shows. 

Finally, we notice that the Go Wo [Gdft] DOS is signifi- 
cantly more broadened than both the GoWo[Ghf] and GW 
DOS. The reason for this is that the DFT DOS has a relatively 
large weight close to Ep. This enhances the QP scattering and 
leads to shorter life-times of the QP in the Go Wo [Gdft] calcu- 
lation. Noticing that the QP life-time is inversely proportional 
to Iml^GW this explains the broadening of the spectrum. 



Vni. CONCLUSIONS 

With the aim of investigating the role of electronic corre- 
lations in quantum transport, we have implemented the non- 
equilibrium GW approximation to the electronic self-energy 
of a finite region of interacting electrons coupled to non- 
interacting leads. We have shown, both analytically and by 
means of numerical examples, that the self-consistent GW 
self-energy leads to identical currents at the left and right in- 
terfaces of the central region. In contrast, the widely used 
Go Wo self-energy does not conserve particle number and thus 
violates the continuity equation. More generally we have 
shown that any ^-derivable self-energy will yield identical 
left- and right- currents independent of the basis set applied. 

Using a WF basis we have introduced an effective electron- 



electron interaction which resembles the real space represen- 
tation but is spin-dependent and self-interaction free in the 
WF basis. In general this provides a means for reducing self- 
interaction errors in diagrammatic approaches like the GW 
method. 

The GW method was applied to the Anderson impurity 
model. In equilibrium and T = we found that the self- 
consistent GW approximation describes the width of the 
Kondo resonance well for intermediate interaction strengths, 
U = A and F > 0.5. On the other hand the sidebands of 
the spectral function are always missed in GW. We presented 
non-equilibrium IV curves and discussed the important effect 
of quasi-particle scattering under finite bias which reduce the 
QP life-times leading to a broadening of spectral features and 
significant suppression of the finite bias conductance. Finally, 
we demonstrated that the Go Wo approach can produce severe 
errors including violation of charge conservation and negative 
differential conductance. The errors become more significant 
at higher bias and close to magnetic transition points. 

We investigated the properties of a molecular junction con- 
sisting of a benzene molecule sandwiched between featureless 
leads. To describe the benzene we used a minimal Wannier 
function basis set which was shown to reproduce the exact 
Hartree and exchange matrix elements to within 5%. The cal- 
culated ionization potential in GW was found to be in good 
agreement with the experimental value. A significant reduc- 
tion of the GW HOMO-LUMO gap was observed for increas- 
ing molecule-lead coupling. The effect comes from the corre- 
lation part of the GW self-energy and reflects the more effi- 
cient screening in a strongly, compared to a weakly, coupled 
junction. 

Finally, the non-equilibrium differential conductance of the 
benzene junction was calculated in DFT-PBE, HF, and GW 
as well as GoWo[HF] and Go Wo [DFT]. It was found that HF 
and GoWo[HF] yield results similar to GW, while both DFT 
and Go Wo [DFT] yield significantly larger conductances. In 
particular, this shows that the Go-dependence of the GoWq 
approximation should not be disregarded. The trends in con- 
ductance were explained in terms of the size of the HOMO- 
LUMO gap of the molecule which also shows significant vari- 
ation depending on the approximation used. 
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APPENDIX A: HARTREE AND EXCHANGE POTENTIALS 

In this work the exchange and Hartree self-energies have 
been evaluated from Eqs. i45[ and (l46T l with the Coulomb 
matrix elements restricted to a certain subset (the set V^^'^ 
defined below). Here we investigate the quality of such ap- 
proximations by testing their ability to reproduce Hartree 
and exchange energies of the molecular orbitals of a ben- 
zene molecule^l. We thus consider the following truncation 
schemes 

V^'^ = (Al) 

V^^^ - V[{V,J,,J},{V,J,,,}AVu,JJ},{V^^.^J}i (A2) 

1/(3) = v[ {v,j,,,}, {Vu^jj}, {v,,,^J}, {v,k^jk(m) 

where e.g. the notation ^^[{Vij.y }] means that all elements of 
the form Vij^ij are included in the sum in Eq. ([3]l. 

The molecular orbitals of benzene, {V'n}. can, by construc- 
tion of the WFs {<j)i}, be exactly expanded as, 

V'nW = ^Ci„(/),(r). (A4) 

i 

I 



The 1 8 WFs used to describe the benzene molecule are plot- 
ted in Fig. [TTtb). For the molecular orbital tpn we can then 
calculate the exact Hartree and exchange energies from 

= 2f:/drdr' ^"(^)^^'"f^^^*t7^'^^''"^'^ 
Alternatively we can insert the expansion (IA4b and get 

{'ipn\^h\lpn) = '^C*^T,h,ijCjn (A5) 

ij 

{i>„\^x\i>n) = ^C*„T,x^ijCj„ (A6) 
ij 

where J^h.ij 'Sx,ij are the self-energies in the WF basis 
obtained from Eqs. (I45]l.(l46b. The latter are approximated by 
the truncation schemes ( IA1HA3I ) for the Coulomb integrals, 

Vij,kl- 

In table U we compare the exact values of the Hartree and 
exchange matrix elements for the frontier molecular orbitals 
to the approximate ones obtained using the truncated interac- 
tions. We note that V'^'^\ which is the truncation scheme we 
have used, leads to average deviations around 5%. 

As a final remark we notice that our results for |Ea; | ■(/'„) 
evaluated using V"(^) provides roughly the same accuracy as a 
recently developed method combining tight-binding DFT with 





edft 




(tpnl'Sxl'tpn) 


State (symmetry) 






1/(2) 


1/(3) 


Exact 




1/(2) 


1/(3) 


Exact 


HOMO-2 (tt) 


-8.94 


217.6 


221.2 


233.2 


233.0 


-14.6 


-16.4 


-16.3 


-16.7 


HOMO-1 (a) 


-8.12 


253.6 


244.7 


224.7 


224.6 


-24.7 


-20.2 


-20.2 


-19.6 


HOMO (tt) 


-6.20 


220.1 


223.0 


229.5 


229.5 


-13.7 


-15.2 


-15.0 


-15.1 


LUMO (tt*) 


-1.08 


222.1 


222.7 


219.1 


219.3 


-7.2 


-7.5 


-7.5 


-7.2 


LUMO+1 (tt*) 


2.68 


223.1 


221.4 


199.3 


199.7 


-6.1 


-5.3 


-5.8 


-4.7 


Average deviation (%) 


7.3 


5.8 


0.1 




15.5 


4.5 


6.7 





TABLE I: Hartree and exchange energies in eV for five frontier molecular orbitals of the benzene molecule. The values are obtained using 
the truncated interactions defined in Eqs. ( IAHA3t as well as the full interaction V (the exact result). For reference the first column shows the 
eigenvalues as calculated using the PBE xc-functional. 



APPENDIX B: ASSESSMENT OF EFFECTIVE 
INTERACTION 

As discussed in Sec. lIII Al the GW approximation includes 
only a single diagram at each order of the interaction. The 
error resulting from such an approximation is - to lowest order 



- similar to the error of approximating a HF calculation by a 
Hartree calculation. It is not obvious that the best result of 
such an approximation is obtained by using the full interaction 
of Eq. (|3]l. For example such a strategy would lead to self- 
interaction errors. 

In table (middle panel) we compare Hartree matrix el- 



21 



ements of some molecular orbitals of benzene^ , evaluated 
using different effective interactions. Notice that the values 
listed in the two leftmost columns differ by the inclusion of 
the spin dependent term of Eq. dSTT l. In the right column we 
show the exact HF result, i.e. the correct result to first order 
in the interaction. The last row shows the average deviation of 
the Hartree energies from the exact HF energies. 

From table|II]we conclude that the effective interaction pro- 
duces results of comparable accuracy to the full interaction, if 
one attempts to reproduce the exact result to first order from 
the Hartree approximation only. The fact that T4if performs 
better than {Vij,ij} indicate that the spin-dependent term in 



t^ff, which removes the self-interaction in the WF basis, is 
significant. 

Extrapolating these observations to higher order we con- 
clude that the use of V^s in GW calculations should produce 
results comparable to GW calculations based on the full in- 
teraction. 

At this point we stress again, that for practical calcula- 
tions we use the truncation scheme of Eq. ( IA2l i for evaluat- 
ing Hartree and exchange. Thus the results presented in this 
section only serve to estimate the performance of the effective 
interaction for the higher-order GW diagrams. 







(V'n|Eh + E^|V'n> 


State (symmetry) 


{Vij.ij} 




Exact 


Exact 


HOMO-2 (tt) 


217.6 


207.4 


233.0 


216.3 


HOMO-1 (a) 


253.9 


230.1 


224.6 


205.0 


HOMO (tt) 


220.7 


210.1 


229.5 


214.4 


LUMO (tt*) 


222.8 


212.2 


219.3 


212.1 


LUMO+1 (tt*) 


223.7 


213.1 


199.7 


195.0 


Average deviation (%) from exact HF 


9.5 


5.5 


6.0 




(right column) 











TABLE II: Left part: Hartree self-energy for some of the frontier orbitals of the benzene molecule. The Hartree self-energy has been evaluated 
using the effective interaction Eq. OOt . the effective interaction without the spin dependent correction (second term in Eq. ll31t), and using the 
full interaction Eq. ([3j (exact result). Right: The exact value of the Hartree-Fock self-energy. Note that the spin dependent correction term in 
Vcff cancels the self-interaction (in the local Wannier basis) and thus incorporates part of the exchange in the Hartree potential. Last row shows 
the average deviation of the Hartree potential from the exact Hartree-Fock potential. 



APPENDIX C: A USEFUL RELATION 

Let B{t, t') and C(r, r') be two matrix valued functions on 
the Keldysh contour, and consider the commutator A defined 
by 

A{t,t') = y[S(r,ri)C(Ti,r') -C(r,ri)B(Ti,r')]dn, 

(CI) 

where matrix multiplication is implied. Under steady state 
conditions where the real time components of B and C can 
be assumed to depend only on the time difference t' — t, the 
following identity holds: 

Ti[A<{t,t)]^ J ^Tr[i3<HC>H-i?>HC<H]. 

(C2) 

To prove this relation we first use the Langreth rules to obtain 



- C<{t,h)B''{h,t') - C''{t,h)B<{h,t')]dti. 



Since all quantities on the right hand side depend only on 
the time difference we identify the integrals as convolutions 
which in turn become products when Fourier transformed. We 
thus have 



2^ 



[B<{lu)C''{uj) + B'-{uj)C<{lo) 



-C<{lu)B''{lu) - C''{uj)B<{uj)]. 

Eq. iCH now follows from the cyclic property of the trace and 
the identity G'' - G" ^ G> - G< . 



APPENDIX D: COUPLING TO QUASI-PERIODIC LEADS 

We consider the coupling of the central region (C) to the 
left lead (i) in the case where L is periodic only beyond a cer- 
tain transition region (T). We refer to the periodic parts of the 
lead as principal layers and denote the corresponding blocks 
of the Hamiltonian matrix by Hq. Without loss of generality 
we assume nearest neighbor coupling between the principal 



22 



layers and denote the coupling matrices by vq. The transition 
region is assumed so large that there is no coupling across it, 
i.e. between the central region and the first principal layer If 
this is not the case the transition region must be extended by 
the first principal layer The Hamiltonian of the left lead and 
its periodic part can then be written as 



V 



ho 




vo 
ho 



: \ 



vt 
hr ) 



I 



V 



'^o 

t 

■"0 



: \ 



ho ) 





(Dl) 

The (retarded) GFs defined from and hp^^ are denoted by 
50, L and .gg"^^, respectively. The lower right block of (?o,Lj cor- 
responding to the transition region, is denoted [.go.Llr, and the 
lower right block of g^"^, corresponding to the first principal 
layer, is denoted by [Sq'lIo- ^^^^ '■^^ following equation 

[ffO,L]T = [(^ + «7?)/-/jT-ST]"' (D2) 

where the self-energy is given by 

St = Wrl^o^Lln'^r- (D3) 

In the above equation [s'q'^^Jo can be obtained using the stan- 
dard decimation technique^^. The coupling self-energy 
can now be constructed from [go, l]t and the matrices Htc 
and hcT which describe the coupling between the transition 
region in the left lead and the central region. 



= hcT[9Q._ 



TC- 



(D4) 



We remark that her and hxc are sub-matrices of hcL and 
h^c- Completely analogue results hold for the coupling to 
the right lead. 



^0 



Central region 



h. 



h. 



Periodic part 



Transition region (T) 



Left lead 



FIG. 16: The transition region, T, is defined as the part of the lead 
beyond which the lead Hamiltonian becomes periodic. 



APPENDIX E: RETARDED FUNCTIONS FROM 
CORRELATION FUNCTIONS 



In steady-state all four real time GFs Eqs.(fT6ll-(fT9]l follow 
from the retarded and lesser components and thus it suffices 
to calculate these. 

Given G'"(w) and G^((jj) sampled on an equidis- 
tant frequency grid, the corresponding GW self-energy, 
SGw[G](ti;), can be obtained from Eqs. (l36]l-(l4Tli using the 
fast Fourier transform (FFT) to switch between energy and 
time domains. However, as alternative to Eqs. (l36l ) and (l40l) 
we have found it more useful to obtain Sqj^? and from the 
relation 



x^{t)^e{-t)[x>{t)-x<{t)i 



(El) 



which is valid for any function X on the Keldysh contour that 
does not contain delta functions. Note, that when applied to 
TiQW Eq. ( lElb yields only the correlation part of Sg^^- as ex- 
plained in Sec. nil B1 The reason why (lEll) is so useful is that 
X^ [lo) falls off as 1 /lo (due to the step function in time) which 
makes it difficult to obtain a faithfull representation of X"^ [t) 
from an FFT of X'^{ljj). In contrast are well lo- 

calized (they are smooth in time), and the FFT can be safely 
used to obtain X'^/^{uj) from and vice versa. It is 

possible to reduce the size of the frequency grid significantly 
if a zero-padding of X^/''{uj) is introduced before the FFT 
is applied to obtain As discussed in Sec. IIIIBl 

Eq. IE II with X = H yields the correlation part of the GW 
self-energy. The static Hartree and exchange terms, Yih and 
Sj,, are calculated from Eqs. ( |45] l and ( |46] l. Once the self- 
energies have been calculated a new set of GFs can be calcu- 
lated from Eqs. dlHl and (|26] |. 
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